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In terminal airspace, integrating arrivals and departures with shared waypoints provides 
the potential of improving operational efficiency by allowing direct routes when possible. 
Incorporating stochastic evaluation as a post-analysis process of deterministic optimiza- 
tion, and imposing a safety buffer in deterministic optimization, are two ways to learn and 
alleviate the impact of uncertainty and to avoid unexpected outcomes. This work presents 
a third and direct way to take uncertainty into consideration during the optimization. The 
impact of uncertainty was incorporated into cost evaluations when searching for the opti- 
mal solutions. The controller intervention count was computed using a heuristic model and 
served as another stochastic cost besides total delay. Costs under uncertainty were evalu- 
ated using Monte Carlo simulations. The Pareto fronts that contain a set of solutions were 
identified and the trade-off between delays and controller intervention count was shown. 
Solutions that shared similar delays but had different intervention counts were investi- 
gated. The results showed that optimization under uncertainty could identify compromise 
solutions on Pareto fonts, which is better than deterministic optimization with extra safety 
buffers. It helps decision-makers reduce controller intervention while achieving low delays. 

I. Introduction 

In the National Airspace System (NAS), terminal airspace is often busy and complicated because hundreds 
of flights must fly through a limited airspace in a short time period and most flights in terminal areas are 
climbing and descending with varied speeds. Situations can become severe when there are several neighboring 
busy airports. Bottlenecks can be formed easily and impair the efficiency of air traffic operations. Therefore, 
improving the operation efficiency in terminal areas with narrow airspace and shared resources is critical for 
building an efficient air traffic system. 

In past years, to improve the operational efficiency in terminal airspace, some researchers attempted to 
solve arrival scheduling problems efficiently . 1-6 Still others 7-9 addressed airport surface management prob- 
lems, which are correlated to the efficiency of terminal airspace operations. When different arrival and/or 
departure flows share the same resources such as runways, waypoints and/or route segments, inefficient 
operations emerge because of the constraints of shared resources. Such interactions can happen among de- 
partures, arrivals, or between departures and arrivals, which is common in terminal areas. These interactions 
are expected to happen more often when continuous descent approaches are deployed in the future. Recent 
studies 10-13 showed that optimized integrated arrivals and/or departures in major airports or metroplex 
areas have promise for improving operation efficiency. However, the benefits from above optimal schedules 
are calculated under deterministic scenarios and they are usually sensitive to flight time uncertainties, which 
can be caused by many sources, such as inaccurate wind prediction, error in aircraft dynamics, or human fac- 
tors. Therefore, the impact of uncertainty must be taken into account when evaluating a schedule’s benefits. 
Incorporating stochastic evaluation as a post-analysis process of deterministic optimization as in previous 
work 14 is one passive way to learn the impact of uncertainty and to avoid unexpected results. Imposing an 
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extra safety buffer upon separation requirements is another prevalent yet passive way to consider uncertainty. 
Optimizing schedules under uncertainty is a proactive method that directly takes uncertainty into account. 
It provides a set of compromise solutions on the Pareto front and helps decision makers avoid unexpected 
effects from uncertainty. 

This work presents an optimization-under-uncertainty method for integrated arrival and departure schedul- 
ing problems. The impact of uncertainty was evaluated directly in the optimization. The optimization is 
multi-objective including total delay and controller intervention count. Both costs were evaluated stochasti- 
cally using Monte-Carlo simulations. To enable the time-consuming optimization, the number of Monte-Carlo 
simulations were decreased without sacrificing accuracy. The Pareto fronts that contain sets of solutions were 
presented to show the trade-off between delays and interventions. Solutions with similar delays but different 
controller intervention count were investigated. Furthermore, the solutions provided via deterministic opti- 
mization with extra buffers were presented and compared against the Pareto front provided by stochastic 
optimization or optimization under uncertainty. 

In the paper, Section II revisits the problem studied in previous work. Section III presents methodol- 
ogy that includes modeling, optimization objectives, uncertainty evaluation approach and implementation. 
Section III provides the results and analysis. Section IV provides conclusions for this work. 

II. Background and Problem 

The interactions between arrivals and departures in Los Angeles terminal airspace were presented in previ- 
ous work 12 for studying optimal integrated operations in deterministic circumstances. The uncertainty anal- 
ysis of the deterministic solutions has been conducted in other work 14 as determined from post-processing. 
This work is focused on the same problem of integrating arrivals and departures in Los Angeles. 



Figure 1. Interactions between SADDE arrivals and CASTA departures 


According to the Standard Terminal Arrival Routes (STARs) and the Standard Instrument Departures 
(SIDs) of Los Angeles terminal airspace, the arrivals from the FIM fix would follow procedure SADDE6 ( 
FIM-SYMON-SADDE-SMO) and the departures to the North need to follow procedure CASTA2 (RWY- 
NAANC-GHART-SILEX) (see Fig. 1). The arrivals are requested to maintain their flight altitudes above 
10,000 feet at fix GHART and the departures have to keep theirs at or below 9,000 feet at the same fix 
to procedurally avoid potential conflicts between arrivals and departures. If there were no interactions, 
departures to the North and arrivals from FIM would have flown direct routes. As shown in the Fig. 1, the 
direct routes would be RWY-WPT2-WPT1 and FIM-WPT1-SMO for departures and arrivals, respectively, 
where WPT1 and WPT2 are made-up fix names for simplicity. Compared to these direct routes, besides 
flying non-preferred altitudes, individual arrival and departure flights following current procedures will waste 
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approximately 60 and 120 seconds, respectively. 

Table 1 shows a representative schedule of 14 flights, which covers half an hour of traffic in actual 
operations. Two flows are included: 6 departures to the North from Runway 24L (RWY) and 8 arrivals 
from FIM. This schedule used the traffic between 9:00am and 9:30am (local time) on December 4, 2012 as 
reference. The initial times shown in the table are relative times to simulation start time. The “Order” of 
each flight is sorted based on intial times. 


Table 1. Scheduled initial times 


Order 

FIM (sec) 

RWY (sec) 

1 

39 

68 

2 

446 

165 

3 

728 

363 

4 

1106 

529 

5 

1332 

1613 

6 

1475 

1830 

7 

1613 

NA 

8 

1770 

NA 


III. Method 

In this study, the model is constructed similarly to previous works. 12,14 However, the optimization is 
formulated as a multiple objective optimization in order to include more uncertainty-related costs. The 
first objective is to minimize delay and the second objective is to minimize controller intervention, which 
is believed to correspond with controller workload increase when uncertainty arises. Both objectives were 
evaluated using Monte Carlo simulations. 

A. Model 

As the route structure adopted in this study is the same as in previous work, the model that is built upon 
the route structure is also the same. Detailed description of the model should be found in previous work. 12 
For completeness, this model is briefly revisited. Three separation methods were compared in previous 
work: spatial, temporal, and hybrid. Spatial separation uses the same strategy as in SIDs and STARs to 
spatially separate interacting departure and/or arrival flows. Temporal separation utilizes the direct routes 
with conflicts resolved merely with temporal control. Hybrid separation applies both temporal and spatial 
separations. This work was focused on hybrid separation. 

In the formulation of hybrid separation, four design variables were defined for each FIM arrival i: dl(piM,i) 
is the delay before or at FIM; r^ FIM ^ is the route option, where 0 denotes the direct route and 1 denotes 
indirect route; is the aircraft speed between FIM and WPT1 for the direct route or the speed between 

FIM and WPT2 if the indirect route is chosen; d2(F/M,i) is the delay between WPT1/WPT2 and SUTIE 
to ensure separation at SUTIE. For a departure flight j, three decision variables were defined: d^EP,j) is 
the delay before departure; r^EPj) is the route option, where 0 denotes the direct route and 1 denotes the 
indirect route, v^eepj) is the speed from departure to WPT1. 

Separation requirements were applied as hard constraints at fixes that could have potential violations, 
such as FIM, RWY, WPT1, WPT2, and SUTIE. Separation requirements were 3 nmi at the RWY and 4 nmi 
elsewhere. Additionally, an extra safety buffer <5 could be added in the separation constraints if required. 
For example, Eqns. 1, 2, and 3 show separation constraints for crossing flights between FIM arrivals and 
departures at WPT1, where tpjM (wpti,i) is the arrival time of FIM arrivals at fix WPT1 and tp>EP(wPTi.j) 
is the arrival time of departures at fix WPT1. Eqn. 1 showed that if both arrival flight i and departure flight 
j took direct routes the separation must be satisfied. The separation requirement is 4.0 nmi in distance and 
in the time scale the separation depends on the speeds of both flights. 
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(1 - r( F iM,i))-{ 1 - f(DEP,j))-[tFIM(WPTl,i) ~ tDEP (W PTl,j)} 



( 1 ) 


V(FIM,i), if tFIM(WPTl,i) < tDEP(WPTl,j) 
V(depj), otherwise 


( 2 ) 


a = 


atan[ 

atan[ 


V(DEP,j) 

V(FIM,i) 

V(FIM,i) 

V(DEP,j) 


if tFIM(WPTl,i) < tDEP(WPTlJ) 
otherwise 


(3) 


B. Evaluation of uncertainty 

To explain how to evaluate costs caused by uncertainty, three parts are discussed: perturbation error sources, 
heuristic controller behavior model, and identification of uncertainty cost evaluation method. 

1. Uncertainty sources 

There exist many sources of uncertainty: inaccurate wind prediction, unknown aircraft weight, error in 
aircraft dynamics, or inaccurate prediction of flight take-off or arrival time. Because these errors are mainly 
reflected in flight arrival or departure times in this study, to simulate the uncertainty, errors were added in 
flight arrival times at waypoint FIM and SUTIE, and departure times at RWY, respectively. The imposed 
errors were assumed to follow normal distributions. Although, other distribution may be possible, such as 
Poisson in Mueller’s work, 15 it should be straightforward to change in the future when the distribution is 
finally identified. In this work, the arrival time error follows a standard deviation of 30 seconds with a zero 
mean, which was commonly used as a desired and reliable prediction accuracy in arrival trajectory prediction 
studies. 16, 17 The departure time error’s standard deviation is 90 seconds and the mean value is 30 second, 
which is based on the Call For Release (CFR) 3- minute compliance window. 18 The CFR window is often 
structured to allow departure 2 minutes prior to or 1 minute later than the target coordinated departure 
time. 

2. Heuristic controller behavior model 

A heuristic model was built to mimic controller intervention behaviors as in previous works. 12,14 When 
stochastic errors are added in flights’ entry times, the heuristic model resolves potential conflicts by imposing 
extra delays to corresponding aircraft, while keeping the same route options as in the given solution. The 
conflict resolution simulates controller’s behavior by following the First-Come-First-Served (FCFS) rule. At 
a waypoint, this model sorts the sequence based on the flights’ perturbed entry times. Then extra delays are 
imposed if the separation between any two adjacent aircraft doesn’t meet the defined requirement, which is 
the sum of the separation requirement and an additional buffer if required. The imposed delays are then 
propagated to following waypoints. In actual operations controllers might use improved heuristics other than 
this strategy, it should not be a problem to incorporate it into this formulation in the future. 

3. Evaluation method 

A traditional approach to evaluate a function with random inputs is the Monte Carlo method, which gen- 
erates random realizations for the prescribed random inputs and utilizes repetitive deterministic function 
evaluations for each realization. In this study, the deterministic function is the heuristic controller behavior 
model. A brute-force Monte Carlo method with K realizations converges asymptotically at a rate 1 /y/K. 19 
Groups of methods in the literature have been proposed to speed up the evaluation process while keeping 
similar accuracy. Many methods were developed for solving stochastic differential equations with uncertain 
inputs, such as generalized polynomial chaos (gPC), 20 multi-element generalized polynomial chaos (ME- 
gPC), 21 and multi-element probabilistic collocation method (ME-PCM). 22,23 The probabilistic collocation 
method (PCM) has been applied successfully to dynamic simulation of power system 24 and geophysical 
models. 25 It has also been tried in certain air traffic problems. 26,27 The main idea of these methods is 
to speed up the process by applying a smart sampling method with small samples. This is often achieved 
by a certain type of decomposition which can approximate the random process with desire accuracy. For 
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instance, PCM uses a polynomial expression to represent the random response as a function of orthogonal 
polynomials of random variables or polynomial chaos expansion. However, although these methods are good 
at solving stochastic differential equations or time-step simulations, limitation exists when the dimension of 
the random inputs is high and/or regularity of the function is poor. Experiments were conducted for the 
problem in this study using PCM, ME-gpc, and MEPCM-A methods. But they were not successful. Poor 
response accuracies were obtained for the first-order moments or mean values, and even worse precisions were 
produced for the second-order moments or standard deviations. The reason may be that the problem in this 
work involves high dimension (number of aircraft) and the response is discontinuous with low regularity. 

30 | | | | | | | | | | | | | | 
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Figure 2. Percentage difference for mean values with varied number of Monte Carlo simulations 


Given the fact that the Monte Carlo method is independent of the dimensionality of the random space, 
experiments were also conducted using the Monte Carlo method. Directly incorporating full scale Monte 
Carlo simulations (say 5,000 or 10,000 simulations) into objective evaluation would be computationally 
expensive and make the optimization impractical. However, experiments showed that the calculations of 
mean values and standard deviations converge fast when the number of Monte Carlo simulations goes beyond 
1,000. Figure 2 shows the differences in the controller intervention costs from 10 to 8,000 simulations during 
optimization process for around 10,000 different scenarios/cases. The horizontal axis is the number of 
simulations and the vertical axis denotes the relative differences of mean values against the “truth” , which is 
calculated using 12,000 simulations. The differences are expressed in percentage. It can be seen that at 800 
simulations the differences are reduced within 3%, and the progress stabilizes as the number of simulations 
increase further. Figure 3 presents the changing rates of the differences, or the derivatives of the curve in 
Figure 2. Apparently, the reduction of the differences is slowed down dramatically after 500. Experiments 
with different scenarios/cases also showed similar trends. Although, this don’t serve as an official proof, it 
does show that 1,000 Monte Carlo simulations could provide a good approximate to the random process for 
this work. 


C. Objectives 

The optimization in this work is formulated as multiple objective: minimizing expectation values of total 
delays and controller interventions, respectively. The objectives are shown in Eqn. 4: 


J\ — 'unimpeded 

Ji= EEiJVi] 


(4) 


5 of 11 


American Institute of Aeronautics and Astronautics 


0.7 


0.6 - 


0.5- 


0.4- 


0.3- 


£ 0 . 2 - 


0.1 


°10 


30 50 70 100 500 

Number of simulations 


-V. A i J 1 1 

1000 15002000 3000 5000 8000 


Figure 3. Changing rate of percentage difference with varied number of Monte Carlo simulations 


Ji is the expectation value of the total delay over the sampling space caused by uncertainty. The total 
delay can be expressed as the difference between the sum of flight transit times U and the sum of the 
unimpeded transit times. The unimpeded transit time is the transit time when a flight takes a direct route 
without any delay. As the sum of unimpeded transit times is a constant, its expected value is represented by 
Tunimpeded in the equation. J 2 is the expectation value of the sum of controller intervention count. When a 
flight needs to be delayed due to a separation violation at a merging or diverging waypoint, the intervention 
count will be increased be one. Both costs are evaluated over sampling points using aforementioned Monte 
Carlo simulations. 

The optimization is solved using a Non-dominated Sorting Genetic Algorithm (NSGA) because of NSGA’s 
ability of handling multiple objective optimization . 12,28 As the costs are handled independently, NSGA will 
lead its search to a Pareto front, where no solution on the front is dominated by another. The dominance 
with two objectives is defined as: Assuming the objective is to minimize costs and the constraint function g 
has to be nonnegative, solution A is dominated by solution B only if: 

Ji(A) > Ji(-B), if g A > 0 and g B > 0, or, g A = g B 

< J 2 (A) > J 2 (B), if g A > 0 and g B > 0, or, g A = g B (5) 

< gB , if g A < 0 and g B < 0, or, g B > 0 and g A < 0 

where J\ and J 2 are the objectives and g is the constraint value. 

D. Implementation 

The optimization was initially implemented using C and Monte Carlo was multi-threaded with pthread. 
Using 1,000 simulations, this optimization took about 6 hours to solve the 14 aircraft scenario on a Ma- 
cOS platform with 2x2.66 GHz 6-Core Intel Xeon and 8 GB RAM. The running time is prohibitively long, 

which makes applications or experiments impossible. To improve the performance, in recent work, 29 the 

Monte Carlo part was implemented using CUDA programming with GPUs, and the running-time was re- 
duced to around 2.5 minutes on a Linux platform with one GeForce GTX690 GPU. After further code 
optimization/modification, solving the same scenario now takes around 30 seconds on a Linux platform with 
18x2.5 GHz Xeon and 32 GB memory and two GeForce GTX690. This fast-time implementation makes the 
stochastic optimization a promising method in application. 
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IV. Results 


This section presents the solutions optimized under uncertainty. The trade-offs between delay reduction 
and controller intervention count increase are shown. The solutions that have similar delays and different 
intervention counts are compared. The advantages of stochastic optimization over deterministic optimization 
are shown through results. 

A. The Pareto front 

Given two objectives of reducing delays and controller intervention counts, NSGA led the search to a Pareto 
front, where no solution on the front was dominated by other solutions. GA-based optimizations are sensitive 
to initial guesses, which are decided by initial randomized seeds. Multiple runs are necessary to increase the 
chance of getting optimal solutions. Ten runs were performed in this work. Solutions in final generations of 
the evolution process were recorded as they were close to the Pareto front. Figure 4 presents the solutions 
in final generations ( about 700 solutions) that are on the Pareto front. Each dot (Ji, J 2 ) corresponds 
to a solution and its coordinates J\ and J 2 are the two costs of the solution. The green curve presents 
coordinates/costs calculated in the optimization process, using 1,000 Monte Carlo simulations. Whereas, 
the black curve presents coordinates computed in a post-process using 5,000 Monte Carlo simulations. The 
percentage differences are 1.4% and 0.2% for intervention count and delay, respectively. The marginal 
difference further demonstrates that reducing the number of Monte Carlo to 1,000 works well. The big 
dot at coordinate (1253, 1.23) sets a reference point for intervention counts and delay costs under spatial 
separation, which can be treated as current ATC procedure. To retrieve the best solution under spatial 
separation, an optimization was conducted under deterministic scenarios with delay as the single objective. 



Figure 4. Solutions from stochastic optimization with multiple runs 


According to the final solutions, it is noted that a clear tradeoff exists between delays and intervention 
counts. When delay is reduced the chance of controller intervention increases. Variations of intervention 
counts are large for solutions with similar delays, especially when delays are low. A multiple-objective opti- 
mization that incorporates uncertainty can clearly help to find the target delay with minimum intervention 
or visa versa. If a decision support tool can be built upon this method, in terms of the solutions that are 
close to the Pareto front, decision makers can choose a “compromise” solution according to their preferences. 
For instance, if controller intervention is believed to be more important than delay savings, the solutions that 
have a delay of about 900 seconds and intervention count similar to the spatial solution should be picked. 
Or, if two costs are weighted similarly, the compromise solutions around the middle (e.g. delay around 500) 
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can be chosen because of their large delay savings and tolerable controller intervention increase. If controller 
intervention doesn’t cause much workload due to advanced equipage in the future, then a solution with 200 
second delay can be chosen. 


B. Effect of multiple objectives 


To get insight into the effect of the optimization objectives, solutions that shared similar delays but with 
different controller intervention counts are retrieved for examination. For instance, at delays around 430 
seconds, the expectation of intervention count can vary from 2.4 to 3.5. The solution with 3.5 intervention 
count won’t be at the Pareto front. Figures 5 and 6 are the time lines for these two solutions at either end 
of this intervention count range. The small gray boxes are the minimum separation requirements associated 
with flights. As mentioned previously, the separation requirement is a function of aircraft speed and the type 
of potential conflicts (crossing or in-trail). Therefore, the gray boxes have different lengths. Any departure 
following the direct route should go through WPT2, whereas any arrival from the FIM fix flying the direct 
route would pass WPT1. 
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Figure 5. 


Time line for the solution with high controller intervention count 
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Figure 6. Time line for the solution with low controller intervention count 


In terms of costs, these two solutions shared similar delays at around 423 seconds but resulted in different 
controller intervention efforts. The solution associated with Fig. 6 has an expected intervention count of 2.45, 
whereas the Fig. 5 solution corresponds an expected intervention count of 3.56 - a 45 % increase. Based on the 
comparison, the major differences between them are that: The high-controller-intervention solution allowed 
FIM003 and FIM008 to fly direct routes whereas the low-controller-intervention solution chose FIM007 and 
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DEP006 to fly direct routes (See the flights pointed by arrows in the figures). Those differences in route 
options would not affect the delay reduction but would lead to different controller intervention. This subtle 
difference could not be easily predicted without the optimization under uncertainty. 

C. Stochastic vs. deterministic 

When dealing with uncertainty, one prevalent method in scheduling problems is to add an extra buffer on 
top of separation requirement in a deterministic optimization. The choice of the buffer size depends on 
applications and/or users. One advantage of optimization based on this method is its fast running time, 
because it is a deterministic optimization. However, this method doesn’t include any distribution knowledge 
of error sources, simply and blindly enforcing a constant buffer doesn’t usually help to find “best” solutions. 



Delay (s) 


Figure 7. Stochastic optimization versus deterministic optimization 


In order to demonstrate the difference, Fig. 7 presents the comparisons between stochastic and deter- 
ministic optimizations. Because the intervention count cost is zero for a feasible solution in deterministic 
cases, optimization in deterministic scenarios can only contain one objective - the delay cost. The curve on 
the lower left side represents the Pareto front generated using stochastic optimization or optimization under 
uncertainty. Three solutions on the top left side are generated using deterministic optimization with hybrid 
separation and varied buffer size. The circular dot is the solution with Os buffer, and the triangle and square 
solutions are produced with 30s and 60s buffers, respectively. Apparently, the deterministic solutions are not 
at the Pareto front, which means, with the same delay, stochastic optimization can find better solution with 
lower controller intervention cost than deterministic optimization. The solution’s distances to the front are 
random, which implies that the size of the buffer produces a random effect on the solution quality. To be 
complete, Fig. 7 also provides the solutions generated under deterministic optimization and spatial separa- 
tion with varied extra buffers. It is noted that they are overlapped with each other with pretty close statistic 
costs. The reason could be that the natural gaps among flights in this scenario with spatial separation are 
either too big or too small such that imposing extra buffers doesn’t affect the final solution too much. This 
also shows that the effect from the extra safety buffers in deterministic optimization on final solutions can 
not be easily predicted and the usage is blind. As a comparison, with information of error sources, stochastic 
optimization can lead the optimization towards the Pareto front with good-quality solutions. 
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V. Conclusions 


In order to directly take the impact of uncertainty into account, an optimization method of integrated 
arrivals and departures under uncertainty was developed in this study. A simple controller behavior model 
was incorporated to compute controller intervention, such that the controller intervention count could be 
included as another stochastic cost besides delay cost. The problem was then formulated as a multiple objec- 
tive optimization with delays and intervention count costs. Monte Carlo simulations were utilized to evaluate 
stochastic costs, because, unlike other prevalent PCM methods, it is independent of the dimensionality and 
regularity. To enable the time-consuming optimization, the number of Monte Carlo simulations was reduced 
without losing much accuracy. The Pareto fronts that contain sets of solutions were identified. The trade-off 
between delays and controller intervention counts was shown and solutions with similar delays but different 
intervention effort were investigated. Comparison between stochastic and deterministic optimization was 
conducted. As the Monte Carlo simulations and the NSGA algorithm are all well-suited for parallelization 
because of their independent calculations and memories, the current implementation with GPUs can solve 
a typical half an hour scenario in about 30 seconds. 

Through this study, the optimization under uncertainty for integrated arrivals and departures was found 
to be feasible with simplified Monte Carlo simulations. Decreasing the number of simulations from 5,000 to 
1,000 affected the controller intervention evaluation about 3% or less with given samples, but helped reduce 
the computational times to a reasonable level. Using this formulation, the method can provide a sweep of 
solutions at the Pareto front so the decision makers can choose in terms of their preferences. Solutions may 
have similar delays but cause quite different controller intervention efforts. The subtle difference in solu- 
tions/delays can result in a significant difference in intervention counts (e.g. 45%), which would not be easily 
foreseen without the optimization under uncertainty. By the comparison between stochastic optimization 
and deterministic optimization with a safety buffer, it is shown that stochastic optimization outperforms the 
latter by providing better solutions located at the Pareto front. Imposing a safety buffer in a deterministic 
optimization actually produced a random or unpredictable effect on final solutions. 
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